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A NUMERICAL METHOD OF CALCULATING THE BOUNDARY -INDUCED 
INTERFERENCE IN SLOTTED OR PERFORATED WIND TUNNELS 
OF RECTANGULAR CROSS SECTION 

By James D. Keller and Ray H. Wright 
Langley Research Center 

SUMMARY 

A numerical method is presented for calculating the boundary -induced interference 
at subsonic speeds in wind tunnels of rectangular cross section with slotted or perforated 
walls. The slot width or wall porosity can vary throughout the test section. The inter- 
ference can be computed at any point in the test section. The model can be any configura- 
tion and can be located anywhere in the test section. Several examples are given, and 
comparison is made with other methods where available. 

INTRODUCTION 

One of the major problems encountered in the design of a subsonic wind tunnel is 
that of determining a suitable test -section configuration to reduce the interference due 
to the tunnel walls. The use of ventilated (slotted or perforated) wind-tunnel walls has 
proven to be an effective means of accomplishing this goal (ref. 1). Theoretical methods 
are presently available for predicting the interference due to the tunnel walls in certain 
cases. These methods are limited to inf inite -length test sections with constant slot 
width or wall porosity. Some are limited as to model size, position, and load distribution. 

In this report a numerical method is described for theoretically determining the 
boundary -induced interference in ventilated wind tunnels of rectangular cross section. 

The method consists of dividing the tunnel walls into a number of rectangular elements, 
each of which is represented by a source distribution. The boundary conditions are 
satisfied at the centroid of each element. The method considers test sections of finite 
length and can be used to satisfy a variety of boundary conditions. The boundary condi- 
tions need not be constant; thus, varying slot width or wall porosity can be treated with 
this method. 

Several examples are presented, and a sample computer program used in making 
the calculations is given in an appendix. 



SYMBOLS 


a effect of one element on another 

b effect of model on element 

d distance between slot centers 

L lift 

AL/L weighting factor for a lift element 

l slot parameter 

n direction normal to wall 

R porosity restriction factor 

s wing span 

t slot width 

Wj. upwash velocity caused by tunnel walls 

x,y,z Cartesian coordinates 

r m circulation of model 

6 lift interference factor 

£,77,£ Cartesian coordinates 

a source distribution strength 

cp perturbation velocity potential function 

cp* velocity potential function for an element divided by a for the element 


2 



Subscripts: 


i ith element 

j jth element 

m model 

t tunnel walls 


ANALYSIS 


General Statement of the Problem 


The governing equation used in the analysis of the low-speed wind-tunnel interfer- 
ence is 


d%2 9y2 g z 2 


( 1 ) 


where cp is the perturbation velocity potential function for the entire flow field. Let 
cp = cp m + <p^ where cp m is the potential function of the disturbances due to the model 
in free air and cp^ is the potential function of the additional flow due to the tunnel walls. 
If cp m is taken to be a known solution of equation (1) which approximates the flow field 
at a distance from the model in free air, then <p^ can be determined by the fact that cp 
must satisfy certain boundary conditions at the tunnel walls. The objective in deter- 
mining <p^ is to be able to calculate the change in the free -stream conditions caused by 
the tunnel walls. Since cp m needs to be known only on the tunnel walls, any inaccu- 
racies in cf> m near the model will have a negligible effect on the determination of <p^. 


Boundary Conditions 


The boundary condition to be satisfied at a solid wall is that there can be no flow 
through the wall, that is 




where n is the direction normal to the wall. The boundary condition to be satisfied at 
an open jet boundary is that there is no pressure difference across the boundary. This 
boundary condition can be approximated by (ref. 2) 
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Reference 3 shows that the mixed open and closed boundary conditions for a tunnel wall 
with several longitudinal slots can be replaced by a homogeneous boundary condition for 
an ideal slotted wall. This boundary condition is 


cp+ l ^.= 0 


( 2 ) 


where l is a slot parameter given by 

i=2lncsc(l!-') 

7 r \2 d ) 

where t is the slot width and d is the distance between slot centers. Reference 4 
gives the boundary condition for an ideal perforated wall as 


®2 + i. o 

ax R an 


(3) 


where R is an experimentally determined restriction factor. 


Representation of the Tunnel Walls 

In order to satisfy the homogeneous boundary conditions for ventilated (slotted or 
perforated) wind tunnels, the tunnel walls are divided into longitudinal strips and each 
strip is divided into a number of rectangular elements. The boundary conditions are 
satisfied at the centroid of each element. The coordinate system used has the X-axis 
extending along the tunnel center line, with the positive direction being the tunnel stream 
direction. The Z-axis is positive upward, and the Y-axis is chosen so that the coordinate 
system is a right-handed system. Each tunnel wall element is represented by a constant- 
strength source distribution over the element. If cp* is the velocity potential function 
for a particular element divided by the source strength <7 for that element, then 

N 

n - £ n*°j 

i=i 

where N is the total number of elements. 

Consider first an element in the top or bottom wall with corners located as shown 
in figure 1. The potential function at a point (x,y,z) due to this source distribution is 


^2 n d?? dg 

^ ]j(x - |) 2 + (y - V) 2 + (z - ?i ) 2 


(4) 
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Figure 1.- Schematic of an element in top or bottom wall. 

In order to satisfy the boundary conditions, this potential function and its derivatives must 
be evaluated. For convenience in writing the equations, let x^ = (x - x 2 = ( x - £2)? 
Yj * (y - (j ~ ^2)’ and z l * ( z ' ?i)- The required equations are then 
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In order to find the effect of an element in a side wall, it is necessary to interchange 
y and z in equations (5) to (8). 

Computation of Source Strength Slopes 

In order to compute the source strength slopes required to satisfy the boundary con- 
ditions at the centroid of each element, a matrix equation is needed which expresses these 
boundary conditions. Let a^ be the effect at the centroid of the ith element due to the 

source distribution corresponding to the jth element fa = (p* + 1 or a = ^ ~q^~J • 

( 9< ?m 

Let b. be the effect of the model at the centroid of the ith element lb = cp m + 1 ■ — or 
3 (p 1 3 cp \ 

h = HI + 1 EH). Then the matrix equation that expresses the boundary conditions is 

9x R 9n / 

AS = -B ( 9 ) 

where 

A - ^ 

s ■ M 
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B = I b 


and 


Equation (9) can be solved for the cn values which can then be used to compute the 
interference due to the tunnel walls. 

EXAMPLES 

As a relatively simple example, consider the lift interference due to a small lifting 
wing mounted in the center of a square test section. The wing is represented by a horse- 
shoe vortex of circulation T m . The span s of the horseshoe vortex is assumed to be 
so small that it becomes a vortex doublet located at (0,0,0). The perturbation velocity 
potential function cp m due to this representation of the model is given by 


<Pm = 


r m s z 
4 IT y2 + z 2 


<1 + 


Jx 2 + y 2 + z 2 y 


so that 


dcp 


m 


r m s 


47T 


(x 2 + y 2 + z 


2 ) 


3/2 


and 


dcp. 


m 


r m s 9y 


47r(y 2 + z 2 ) 


2yz + 2x 3 yz + 3xy 3 z + 3xyz 3 


(x 2 + y 2 + z 2 ) 


3/2 


9 cp. 


m 


r m s 8z 


47r^y 2 + z 2 ) 


,3.2 


,3_2 


_ 2_2 


y2 _ z 2 + x Y_ + xy* - x°z - xy*z - 2xz 


(x 2 + y 2 + z 2 J 


3/2 


These quantities are used on the right-hand side of equation (9) which is then solved for 
ojT m s. The 

a<p t > 


^y/r m s * The a/r m s values are suitable for the computation of the upwash velocity 


/ w t 

V r m^ 


r m s 9z 


at any point in the test section by summing the velocities due to each 


element. The lift interference factor (ref. 3) is then 


. _ 1 w t 
” 2 r m s 
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Each tunnel wall is divided into four strips of equal width and each strip is divided 
into 10 elements by cutting planes at x = -1.00, -0.70, -0.45, -0.25, -0.10, 0.00, 0.10, 
0.25, 0.45, 0.70, 1.00. (For convenience, the tunnel width is taken to be unity.) In order 
to compare the results with those obtained by using the method of reference 5, let the test 
section have four equally spaced slots in the top and bottom walls and let the side walls be 
solid. Figure 2 shows the lift interference factor at the center of the tunnel as a function 



Ratio of slot width to distance between slots, t/d 

Figure 2.- Lift interference factor for a small-span 
wing in a square tunnel with four slots in the top 
and hottom walls. 

of the ratio of slot width to the distance between slot centers. Also shown are the results 
obtained by using the method of reference 5 for the same case. It can be seen that the 
agreement is excellent. 

Figure 3 shows the longitudinal and lateral distributions of the lift interference fac- 
tor due to the small-span model representation in a test section with - = 0.04. 

d 

Present method 

Method of reference 5 



Longitudinal position^ Lateral position, y 

Figure 3 .- Longitudinal and lateral distributions of lift interference factor for a 
small- span wing in a square tunnel with slotted top and bottom walls. ^ = 0.04. 
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One of the important features of the present method is that it can be used when the 
slot width or tunnel porosity varies. Figure 4 shows the effect on the longitudinal distri- 
bution of the lift interference factor when the slot width is varied linearly in the tunnel 
stream direction. It can be seen from figure 4 that by properly contouring the slot width 


t/d =0.04 

t/d varies from 0 to 0.08 
t/d varies from 0.08 to 0 



Figure 4.- Effect of varying slot width on longitudinal 
distribution of lift interference factor in a square 
tunnel with four slots in the top and bottom walls. 

it is possible to reduce the stream wise variations in the lift interference factor. This 
can be important in testing large models if it is desired to have nearly the same inter- 
ference at the tail as at the wing. 

Figures 5 and 6 show the lift interference factor for a small-span wing mounted in 
the center of a square tunnel with slotted side walls and solid top and bottom walls. The 
results obtained by using the method of reference 6 are also shown in these figures. 



Ratio of slot width to distance between slots,t/d 

Figure 5 .- Lift interference factor for a small-span 
wing in a square tunnel with four slots in each 
side wall. 
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Present method 
Method of reference 6 




Figure 6.- Longitudinal and lateral distributions of lift interference factor for 

a small-span wing in a square tunnel with slotted side walls. — = 0.04. 

d 

For the case of the perforated-wall wind tunnel, consider again the same arrange- 
ment of tunnel elements and the same model representation, but let all four walls be per- 
forated so that the porosity restriction factor R is constant. Figure 7 shows the lift 
interference factor at the center of the tunnel as a function of the restriction factor. Also 
shown are the results from reference 7. 



Figure 7-- Lift interference factor for a small- 
span wing in a square tunnel with perforated 
walls . 
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Side walls of test section 



Figure 8.- Representation of sweptback wing. 

The present method is also applicable to other representations for the model. For 
example, consider the case of a large-span sweptback wing which was presented in ref- 
erences 5 and 6. The wing spans 70 percent of the tunnel width and is swept back 35°. 
The wing is represented by lift elements located at points Pj,P2 , . . Pjq on lines of 
35° sweep as shown in figure 8. The coordinates and the assumed lift distribution for 
each of these points are given in the following table: 


Point 

X 

y 

AL/L 

P 1 

0.0246 

0.0351 

0.1342 

p 2 

.0738 

.1054 

.1334 

p 3 

.1229 

.1756 

.1118 


.1721 

.2458 

.0769 

p 5 

.2212 

.3160 

.0437 

P 6 

.0246 

-.0351 

.1342 

p 7 

.0738 

-.1054 

.1334 

p 8 

.1229 

-.1756 

.1118 

p 9 

.1721 

-.2458 

.0769 

P 10 

.2212 

-.3160 

.0437 


Figure 9 shows the distribution of the lift interference factor along the span of this 
wing for a tunnel with slotted top and bottom walls and also for a tunnel with slotted side 
walls. The same comparison was shown in reference 6 and those results are also shown 
in figure 9. 
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Present method 


Results from reference 6 



Figure 9-- Spanwise distribution of lift interference factor for 

a large-span swept wing in slotted tunnels. — = 0 . 0 6 . 

d 

DISCUSSION 


The present numerical method of calculating wind-tunnel interference has the 
advantage of extreme simplicity. It is simple in concept since it involves only the satis- 
faction of the boundary conditions at discrete points by means of source and sink distribu- 
tions over the elements of the boundary. The mathematics is also simple, involving 
nothing more than algebra, ordinary calculus, and the use of machine computer programs. 
On the other hand, the method has extremely broad applicability, since any model repre- 
sentation is acceptable and the model can be placed anywhere in the test section. A 
further advantage is that the boundary conditions can vary almost without limit over the 
test-section boundaries. This feature should be useful in designing test sections having 
small interference over the whole space occupied by the model, including the wing tips 
and the tail. Another advantage of the present method is the relatively short computer 
time required for a solution. The program presented in the appendix takes only about 
1 minute on a Control Data 6600 computer system. This is in contrast to a run time of 
15 minutes to 2 hours for the programs used in references 5 and 6. It should be pointed 
out that the present method requires only one matrix inversion even though the model is 
represented by several lift elements and the interference is computed at several loca- 
tions in the test section. The methods of references 5 and 6 require that the infinite 
integrals be carried out for each lift element and for each point at which the interference 
is required. The present method could be made even faster in some cases by making 
use of the symmetry of the matrix to be inverted. The present development is oriented 
toward rectangular test sections, but it can easily be adapted to test sections of other 
cross-sectional shapes. 
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The present method has the disadvantage of not permitting the usual assumption of 
infinite length of the test section, but this limitation is regarded as of little importance 
since practical test sections are also limited in length by such considerations as power, 
cost, boundary-layer development, and operating convenience. In fact, the present method 
has the advantage of permitting the investigation of the effects of test-section length and 
of the effects of the upstream contraction region and the divergent diffuser entrance 
region. The accuracy is limited by the size of elements into which the test-section 
boundary is divided, but a judicious selection of size and distribution of elements (smaller 
elements in regions nearer the model) should give satisfactory results with a matrix no 
larger than can be inverted on a large computer. If the matrix is larger than can be 
inverted because of a longer test section or smaller elements, the corresponding simul- 
taneous equations can still be solved by iteration methods, although it is doubtful whether 
the additional labor is justified because of the usual uncertainty regarding the actual 
boundary conditions. 


RESUME 

A numerical method of calculating the boundary-induced interference in slotted or 
perforated wind tunnels of rectangular cross section has been presented. The method 
has broad applicability because it allows for a variation in the boundary conditions on the 
tunnel walls and because the model representation is arbitrary. The method also has 
the advantages of extreme simplicity and short computing time required. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., October 18, 1971. 
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APPENDIX 


SAMPLE FORTRAN PROGRAM 


THIS APPENDIX CONTAINS A SAMPLE FORTRAN PROGRAM FOR COMPUTING THE LIFT 
INTERFERENCE FACTOR IN A WIND TUNNEL OF RECTANGULAR CROSS SECTION WITH 
SLOTTED PERFORATED WALLS* IT IS INTENDED ONLY AS A SAMPLE* MODIFICATIONS 
MUST BE MADE TO THE PROGRAM IN ORDER TO COMPUTE DIFFERENT CASES. THE 
PROGRAM WAS WRITTEN FOR USE ON CDC 6000 SERIES COMPUTERS. 


PROGRAM A3307 C INPUT .OUTPUT .TAPE 5= INPUT »TAPE6=OUTPUT I 

DIMENSION XI (160) »ETA( 160) » ZET A ( 160 ) » X 1 1 ( 160) . XI2( 160) . ETA1( 160 ) . 

* ETA2 ( 163 ) • Z ETA 1C 160 ). ZETA2( 160 . A( 16C. 160) »B( 160) » SIGMA ( 160) 
DIMENSION Cl (160) ,C2( 160 ) ,C3( 160 ) 

DIMENSION XA(10)»YA(10).WT(13) 

SGNl X)=SIGN( 1.0, X) 

P(X1,X2,Y1,Y2,Z> = X2*AL0G(A8S((YIFSQRT(X2*X2«-Y1*Y1*Z*Z>)/(Y2*SQRT( 

* X2*X2+Y2*Y2+Z*Z))) )+Xl*ALOG( ABS ( ( Y2+SQRT( X1*X1+ Y2*Y2+Z* Z ) ) / ( Y1+ 

* SQRT(X1*X1*Y1*Y1*-Z*Z) ) ) ) + Y2*AL0G( ABS( ( Xl+SQRT ( X1*X1+Y2*Y2+Z*Z ) ) t ( 

* X2+SCRT ( X2* X2+Y2*Y2 *1*1 ) ) ) ) + Yl*ALOG( ABS( ( X2+SQRT (X2*X2*Y1*Y14Z*Z) 

* )/( X1+SQRT( X1*X1+Y1*Y1+Z*Z) ) ) > 

* +ABS(Z)*(ATAN(X2*Y2/ABSm/SQRr(X2*X2 + Y2*Y2*-Z*Z) )-ATAN( X1*Y2/ABS( 

* Z)/SQRT( X1*X1+Y2*Y2+Z*Z ) l-ATAN( X2*Y1/ABS(Z) /SQRT(X2*X2+Y1*Y1+Z*Z) 

* )fATAN(Xl*Yl/ABS(Z)/SQRT(Xl*XlfYl*Yl*-Z*Z» ) ) 
S(X1,X2,Y1,Y2,Z)=X2*AL0G(ABS((Y1+SQRTCX2«‘X2+Y1*Y1 + Z*Z))/(Y2+SQRT( 

* X2*X2+Y2*Y2+Z*Z) ) ) ) +X1*AL0G ( ABS ( ( Y2+SGRT ( XI *Xl+Y2*Y2*Z*Z ) ) 1 ( Y1 + 

* SORT ( X1*X1 +Y1* Y1+Z*Z ) ) ) »+Y2*AL0G( ABS( (Xl+SaRT(Xl*Xlf Y2*Y24Z*Z ) )/( 

* X2+SQRT ( X2*X2+ Y2*Y2+Z*Z ) ) ) )+Yl*ALOG( ABS( ( X2+SQRT ( X2*X2+Y1*Y1+Z*Z ) 

* ) MX1 + SQRK X1*X1+Y1*Y1*Z*Z) ) ) ) 

DPDX( X1.X2.Y1.Y2.Z) =ALOG ( AB S( ( Yl + SORT (X2*X2«-Y1*Y1+Z*Z) )MY2+SQRT( 

* X2*X2+Y2*Y2+Z*Z) ) * ( Y2+SQRT ( XI *X1*Y2*Y2+Z*Z > ) / ( Y1 fSQRT ( X1*X1 f Yl*Yi 

* +Z*Z)))) 

DPDY(X1.X2»Y1»Y2.Z) =DPDX( Y1.Y2.X1.X2.Z) 

DPDZ(X1,X2,Y1,Y2,Z)=SGN(Z)*<ATAN(X2*Y2/ABS(Z)/SQRT(X2*X24Y2*Y24Z*Z 

* ) )-ATAN( X1*Y2/ABS( Z ) /SQRT ( X1*X1 + Y2*Y2+ Z*Z) ) — ATAN ( X2*Y1 / ABS( Z ) / 

* SQRT ( X2* X2< - V1*Y1*-Z*Z))*'ATAN(X1*Y1/ABS(Z)/SQRT(X1*X1*-Y1*Y1*Z^Z) )) 

R0( X ) =R00 4DR0*( X— XI 1(1)) 

EL(X)=AL0G(1.0/SIN(R0(X)/2.*PI) )/4./PI 
PI=3. 1415926 


1 

2 

3 

A 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 
29 


THIS PART OF THE PROGRAM DEFINES THE TUNNEL GEOMETRY. HERE IT IS SET UP FOR 
A SQUARE TUNNEL OF UNIT WIDTH ANO HEIGHT. EACH WALL IS DIVIDED INTO FOUR 
STRIPS AND EACH STRIP IS DIVIDED INTO TEN ELEMENTS. THE TUNNEL LENGTH IS 
TWICE ITS WIDTH. 


DO 1 1=1,160,10 
XI1CI)=-1.0 
XI1( 1+1)=-. 7 
XII ( I +2 ) =-.45 
XIl(I+3>=-.25 
XIK 144)=-. 1 
XI1( I45)=0.0 
XI1(I46)=.1 
XI 1( I 47 )=• 25 


30 

31 

32 

33 

34 

35 

36 

37 

38 
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X 1 1( I+8)=.45 

39 


XU(I49) = .7 

40 

1 

CONTINUE 

41 


DO 2 1=1,160 

42 


XI2( I )=XI1( 1+1) 

43 

2 

CONTINUE 

44 


DO 3 1=10,160,10 

45 


XI2( I )=1.0 

46 

3 

CONTINUE 

47 


DO 4 1=1,40 

48 


ETA1I I)=.5 

49 


ETA2C I)=.5 

50 


ETA1 ( I +40) =-.5 

51 


ETA2( 1+40)=-. 5 

52 


ZETAlCI«-80)=.5 

53 


ZET A2 < 1 + 80 ) = • 5 

54 


ZETAK 1+120)=-. 5 

55 


ZETA2I I +120) =-.5 

56 

4 

CONTINUE 

57 


DO 5 1=1,10 

58 


ZETAK I ) = . 25 

59 


ZETA1(I+10)=0.0 

60 


ZETAK 1+20)=-. 25 

61 


ZETAKK30)=-.5 

62 

5 

CONTINUE 

63 


DO 6 1*1,40 

64 


ZETAK K40) = ZETAK I ) 

65 


ETA1I I+80)=ZETA1(I ) 

66 


ETA1I I+120)=ZETA1( I ) 

67 

6 

CONTINUE 

68 


DO 7 1=1,80 

69 


ZETA2II ) = ZET Al( I ) 4. 25 

70 


ET A2 ( 1+80 ) = ET A1 ( 1+80 )+• 25 

71 

7 

CONTINUE 

72 


DO 8 1=1,160 

73 


X I ( I ) = ( XI1( 1 ) +X 1 2 ( I ) ) f 2» 

7* 


ETA(I)=(ETAl(I)+ETA2(I))/2. 

75 


ZETA(I) = ( ZETAK D+ZETA2I I) ) /2. 

76 

8 

CONTINUE 

77 


TSL=XI2I10) 

78 


THIS PART OF THE PROGRAM DEFINES THE WALL CHARACTERISTICS. IN THIS CASE THE 
SIDE WALLS ARE SOLID AND THE TOP AND BOTTOM WALLS EACH HAVE FOUR CONSTANT 
WIDTH SLOTS. THE OPEN RATIO OF THE SLOTTED WALLS IS 6 PERCENT. 


RO0-.06 79 

DR0=0.0 80 

DO 11 1=1,80 81 

C1II)=0.0 82 

C2< I )=0.0 83 

C3 ( I ) = 1.0 8* 

11 CONTINUE 85 

DO 12 1=81,160 86 

C1(I>=1.0 87 

C2 ( I 1 = 0.0 83 

C3(I)=EL(Xim) 89 

12 CONTINUE 90 


15 



APPENDIX 


PRINT 991 91 

PRINT 992 «(I 9 XIl(I)tXI2m«ETAlU)fETA2(I)» ZETA1 (II , ZETA2 ( I ) * 92 

* X I ( I ) » ET A ( I ),ZETA(I),I=1,160) 93 


THIS PART OF THE PROGRAM COMPUTES THE INFLUENCE COEFFICIENTS, A(I,J) 


DO 19C 1=1,160 
NW 1 = 1 

IF(I # GT,4QINWI=2 
IF(I.GT#8C)NM = 3 
IF(I.GT.120)NWI=4 
DO 180 J= 1 , 160 
NWJ = 1 

IF< J.GT*40)NWJ=2 
IF< J.GT.8C)NWJ=3 
IF(J»GT# 120 ) N k J=4 
X1=XI (I)-XUf J) 

X2=X I { I ) — XI 2 f J I 
Y1 = ETA( H-ETAK J) 

Y2=ETA( I)-ETA2< J) 

Z=ZET A ( I)-ZETA(J) 

I F ( NWJ# LT#3)Y1=ZETA(I )-ZETAl(J) 
I F ( NW J# LT« 3 ) Y2= ZE TA (I)-ZET A2 < J ) 
IF(NWJ*LT.3)Z=ETA< I )-ETA(J) 
U=DPDX(X1,X2,Y1,Y2,Z) 
IFCNWUNE.NW J)GO TO 110 
PHI=SCX1,X2,Y1»Y2,Z) 

V=0*0 

I F ( I«EQ*J)V=-2«0*PI 
IF(NWJ.EQ.2)V=-V 
IF(NWJ«EQ*3)W=V 
IF(NWJ*EQ*4> W=-V 
GO TO 130 
110 CONTINUE 

PHI = P (XI, X2,Y1,Y2,Z ) 
V=DPDY(X1,X2,Y1,Y2,Z) 
W=DPDZ(X1,X2,Y1 ,Y2,Z> 

I F < NW J# GT» 2 ) GO TO 120 
T = V 

v=w 

W=T 

120 CONTINUE 
130 CONTINUE 
I F ( NW I • EQ« 

IFCNWI.EO* 

I F (NW I • EQ* 

IF ( NWI • EO* 

180 CONTINUE 
190 CONTINUE 
PRINT 991 
PRINT 993 
* A( I, J+5) , 


94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

I)*V 131 

*V 132 

133 

134 

135 

136 

137 

A ( I, J+3) , A( I , J + 4) , 138 

J+9) , J=1 ,160,10) ,1=1,160) 139 


1 ) A ( I , J )=C1 ( I )*PHI*C2( I )*U*C3( 

2 ) A ( I , J ) =C1 ( I )*PH I +C2( I )*U-C3< I ) 

3) AC I ,J)=C1< I ) *PH I+C2( I )*U+C3m*W 
4 ) A ( I , J ) =C 1 ( I)*PHI+C2( I ) *U-C3 ( I)*W 


,( ( I,J,A( I,J),A(I,J + 1) , A( I , J + 2 ) , 
AC I,J+6),A(l,J+7) ,A( I, J + 8) ,A( I , 


THE FOLLOWING STATEMENT INVERTS THE MATRIX A AND PUTS THE INVERSE IN THE 
PLACE OF THE ORIGINAL MATRIX 

CALL MATRIXC 10, 160, 160,0, A, 160, DET) 140 
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APPENDIX 


THIS PART OF THE PROGRAM COMPUTES THE DISTURBANCE DUE TO THE MODEL, HERE IT 
IS SET UP FOR A NUMBER OF VORTEX DOUBLETS LOCATED IN THE HORIZONTAL CENTER- 
PLANE OF THE TUNNEL. THE PROGRAM FIRST READS THE NUMBER OF LIFT ELEMENTS TO 
BE USED AND THEN READS THE X AND V VALUES AND THE WEIGHTING FACTOR FOR EACH 
ELEMENT. 


READ 994 ,L 141 

READ 995 , ( XA ( 1 1 , YA ( I ) , WT< I ) , I = 1 , L ) 142 

DO 200 1*1,160 143 

BI I )=C.O 144 

200 CONTINUE 145 

DO 299 K= 1 * L 145 

XPP= XA ( K) 147 

YPP=YA ( K ) 148 

DO 210 1*1,40 149 

X=XI(I)-XPP 150 

Y=ET All )-YPP 151 

Z=ZETA ( I ) 152 

PHI=Z/(Y*Y+Z*ZI*< l.C+X/SQRTI X*X+Y*Y+Z*Z)» /4.0/PI 153 

U*Z/4. / PI / ( X* X ► Y*Y «-Z*Z ) **1. 5 154 

V=-(2.*Y*Z+(2.*X**3*Y*Z+3.*X*Y**3*Z*3. *X*Y*Z**3 ) / ( X* XfY*Y fZ*Z > **1. 155 

* 5)/4./PI/(Y*Y+Z*Z)**2 156 

B C I )=BI I)+(C1(I)*PHI+C2( I)*U+C3(I)*V>*WTIK> 157 

210 CONTINUE 158 

DO 220 1*41,80 159 

x=xm»-xpp 160 

Y=ETA(I)-YPP 161 

Z=ZETA(I> 162 

PHI=Z/(Y*Y+Z*Z)*( 1.0+X/SQRTI X*X+Y*Y+Z*ZH /4.0/PI 16 3 

U=Z/4./Pt/<X*XfY*YfZ*Z)**1.5 164 

V=-(2.*Y*Z+<2.*X**3*Y*Z+3.*X*Y**3*Z+3.*X*Y*Z**3I/(X*X+Y*Y+Z*Z)**1. 165 

* 5)/4./PI/(Y*Y+Z*ZJ**2 165 

B(IJ=B(I)HCUI )*PHI fC2I I >*U-C3( I »*V)*WT(K) 167 

220 CONTINUE 168 

DO 230 1=81,120 169 

x=xim-xpp 170 

Y*ETAI I l-YPP 171 

Z*ZET Adi 172 

PHI=Z/( Y*Y+Z*Z> *( 1. 0+X/SQRT ( X*X+Y*Y+Z*Z)» /4.0/PI 173 

U=Z/4./PI/(X*X+Y*Y+Z*Z)**1.5 174 

W= < Y*Y-Z*Z+( X**3*Y*Y+X*Y**4-X**3*Z*Z-X*Y*Y*Z*Z-2.*X*Z**4>/( X*X+Y*Y 175 

* +Z*ZJ**1.5l/4./PI/< Y*Y<-Z*Z»**2 176 

Bm=B(I)+(Cl(I»*PHI+C2l I»*U+C3( I)*W)*WT(K> 177 

230 CONTINUE 178 

DO 240 1*121,160 179 

X*XI(I)-XPP 180 

Y=ETA( I )— YPP 181 

Z = ZETAU) 182 

PHI*Z/(Y*Y+Z*ZI*( 1. O+X/SQRT ( X*X+Y*Y+Z*Z ) > /4.0/PI 183 

U=Z/4./PI/( X*X+Y*Y+Z*Z)**1.5 184 

W=IY*Y-Z*ZM X**3*Y*Y»-X*Y**4-X**3*Z*Z-X*Y*Y*Z*Z-2.*X*Z**4)/( X*X+Y*Y 185 

* +Z*Z)**1.5)/4./PI/(Y*Y+Z*Z)**2 185 

B C I ) =B ( I)+(Cim*PHI+C2m*U-C3m*W)*KT(K) 187 

240 CONTINUE 188 

299 CONTINUE 189 

PRINT 991 190 

PRINT 996 , (E(I I ,1*1,160) 191 
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APPENDIX 


THIS PART OF THE PROGRAM COMPUTES THE SOURCE STRENGTHS WHICH SATISFY THE 
BOUNDARY CONDITIONS. 


DO 3CC 1=1.160 192 

SIGMA(I)=C.O 193 

300 CONTINUE 194 

DO 302 1=1,160 195 

DO 301 J= 1 » 160 196 

SIGMA! I)=SIGMA< I )-A ( I , J ) *B( J ) 197 

301 CONTINUE 198 

302 CONTINUE 199 

PRINT 991 200 

PRINT 996 • (SIGMA ( I ) ,1=1.160) 201 


THIS PART OF THE PROGRAM COMPUTES THE LIFT INTERFERENCE FACTOR. IT READS 
THE X, Y, AND Z VALUES AT WHICH THE INTERFERENCE IS TO BE COMPUTED. 


PRINT 997 202 

400 CONTINUE 203 

DELTA1=0.0 204 

DELTA2=0. 0 205 

DELTA3=0.C 206 

DELT A4=0.0 207 

READ 996 ,XC,YC,ZC 208 

I F { EOF , 5 ) 990 , 40 1 209 

401 CONTINUE 210 

DO 410 J=l,40 211 

X1 = XC-XU(J) 212 

X2=XC-XI2(J) 213 

Y=YC-ETA( J) 214 

Z1 = ZC-ZETA1U) 215 

Z2= ZC-ZET A2 ( J ) 216 

W=DPDY(X1,X2,Z1 ,Z2, Y) 217 

0ELTA1 = DELTAH-W*SIGMAI J>/2. 218 

410 CONTINUE 219 

DO 420 J=41 , 80 220 

X1=XC-XI1(J) 221 

X2=XC-XI2(J ) 222 

Y=YC-ETA( J) 223 

Z1=ZC-ZETA1 C J ) 224 

Z2=ZC-ZETA2( J) 225 

W=DPDY( X1,X2 ,Z1 ,Z2,Y) 226 

OELTA2=DELTA2+W*SIGMA( J)/2. 227 

420 CONTINUE 228 

DO 430 J= 81 » 120 229 

X1=XC-XI1(JI 230 

X2=XC-XI2(J) 231 

Yl= YC-ETA1 ( J ) 232 

Y2=YC-ETA2( J) 233 

Z=ZC-ZETA(J) 234 

W=DPDZ(X1,X2,Y1,Y2, Z) 235 

DELTA3=DELTA3+W*SIGMAC J)/2. 236 

430 CONTINUE 237 

DO 440 J=121,160 238 

X1=XC-XI1(J) 239 

X2=XC-XI2(J) 240 
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APPENDIX 


Y1=YC-ET A1 ( J ) 241 

Y2=YC-ETA2< J ) 242 

Z=ZC-ZETA( J ) 243 

W=DPDZ(X1 ,X2, Y1 ,Y2, Z) 244 

DELTA4=DELTA4+W*SIGMA( JI/2. 245 

440 CONTINUE 246 

DELTA=DELTA1+DE LTA2+DELT A3+DELTA4 247 

PRINT 996 ,XC,YC,ZC. DELTA 248 

GO TO 400 249 

990 STOP 250 

991 FORMAT ( 1H1 ) 251 

992 FORMAT (I4t9F10«6) 252 

993 FORMAT (214*10 FI 0*6) 253 

994 FORMAT (12) 254 

995 FORMAT (3F10«6) 255 

996 FORMAT (10F10. 61 256 

997 FORMAT ( 39H1 X Y Z DELTA) 257 

END 258 
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